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NOMENCLATURE 


constants appearing in equations (23) and (25) 
integral appearing in cross section (see eq. (6)) 
constant associated with intermolecular potential (eq. (23)) 
impact parameter 

constant associated with intermolecular potential (eq. (23)) 
coefficient of wave function (eq. (13)) 
coefficient associated with R 6 term (eq. (23)) 
molecular speed 

/ 3kTT 
~M~ 

transition amplitude 

strength associated with intermolecular potential 
kinetic energy or energy 
distribution function 
relative speed 

Planck f s constant divided by 2 tt 
quantity given in equation (25) 

matrix whose element is defined by equation (20) 

Boltzmann constant 

magnetic quantum number; also molecular weight 

number of simulated molecules within the cell 

number density 

transition probability 

modified transition probability 

collision cross section 



R intermolecular distance 

T temperature 

t time 

V interaction potential 

W eigenenergy of the rotational Hamiltonian 

| a> eigenfunction of the rotational Hamiltonian 
Ax increment of collision time 

6 index of power associated with point centers of repulsion model 

(eq. (5)) or delta function (eq. (16)) 

X quantity defined in equation (25) 

y reduced mass 

v collision frequency 

a effective collision diameter 

t characteristic collision time based on initial translational 

temperature, 1/ (niro 2 ^) 

X deflection angle defined in equation (1) or molecular orientation 

angle (eq. (22)) 

¥ wave function 

a) difference eigenenergy 

Subscripts : 

C distance of closest approach of pair of molecules 

i,j rotational states 

max maximum value 

0 initial 

r rotational 

t total 

tr translational 

a rotational state 
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MONTE CARLO CALCULATIONS OF DIATOMIC MOLECULE GAS FLOWS 


INCLUDING ROTATIONAL MODE EXCITATION 
Kenneth K. Yoshikawa and Yukikazu Itikawa* 
Ames Research Center 


SUMMARY 


The direct simulation Monte Carlo method is used to solve the Boltzmann 
equation for flows of an internally exicted nonequilibrium gas, namely, of 
rotational ly excited homonuclear diatomic nitrogen. In this study, the semi- 
classical transition probability model of Itikawa is investigated for its abil- 
ity to simulate flow fields far from equilibrium. The behavior of diatomic 
nitrogen is examined for several different nonequilibrium initial states that 
are subjected to uniform mean flow without boundary interactions. 

A sample of 1000 model molecules was observed as the gas relaxed to a 
steady state starting from three specified initial states. The initial states 
considered are: (1) complete equilibrium, (2) nonequilibrium equipartition 

(i.e., all rotational energy states are assigned the mean energy level that 
obtains at equilibrium with a Boltzmann distribution at the translational 
temperature), and (3) nonequipartition (i.e., the mean rotational energy is 
different from the equilibrium mean value with respect to the translational 
energy states) . Since only uniform flow is considered, the effect of elastic 
collisions is ignored in the Monte Carlo simulation. 

In all cases investigated the present model satisfactorily simulated the 
principal features of the relaxation effects in nonequilibrium flow of 
diatomic molecules. 


INTRODUCTION 


Understanding the energy balance, as well as the energy transfer 
mechanisms, within the internal excited states of a nonequilibrium rarefied 
flowing gas, is important, in indeed, to a number of problem areas within the 
broad categories of planetary reentry, combustion, and pollution. In this 
paper, results are presented on translation-rotation relaxation. Furthermore, 
these results are obtained by solving Boltzmann T s equation by the Monte Carlo 
direct simulation technique, a method that has received considerable recent 
attention (refs. 1-11). A feature of this method is that it gives insight 
into the effects of relaxation on the microscopic level during molecular 
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collisions; in particular, the instantaneous distributions of internal states 
can be continuously followed. 

The method is described in detail elsewhere (refs. 1 and 6-9). Briefly, 
the flow is computed by following in detail several thousand model molecules 
that are allowed to interact with each other. The coordinates of each mole- 
cule in phase space (including rotational state) are at all times known. 

These coordinates change only during a collision and the modeling of these 
intermolecular encounters is the essence of an accurate simulation. To 
account for these encounters, a molecule and a near neighbor are each selected 
at random as are also their impact parameter and deflection angles - all in a 
manner representative of typical molecules undergoing encounters. They are 
accepted for an interaction or rejected depending on a selected rule that 
depends on cross section and, therefore, on intermolecular potential and the 
relative velocity of the collision pair. Since the initial coordinates (rela- 
tive velocity, impact parameter, and pair of rotational states) are known, 
there remains only to find the final rotational state. This is found by com- 
puting the distribution (transition probability) of all final states accessi- 
ble from the known initial states. The final state is then determined by a 
random selection from this distribution. 

The procedure for handling the n translational n interactions parallels 
other investigations (refs. 1 and 6-9) treating monatomic gases. The proce- 
dure described in this paper differs, however, from other investigations of 
the treatment of internal state interactions. These other investigations fall 
into four categories: (1) semi-empirical, (2) classical, (3) semi-classical, 

and (4) quantum mechanical . The semi-empirical models, energy sink (ref. 11), 
and rough spheres and loaded spheres (ref. 12) , while adequate for steady 
flows at or near equilibrium, lack sufficient physical detail to inspire 
confidence in their use for highly non equilibrium flows. 

The classical models (refs. 13 and 14), although consistent with the 
classical direct simulation Monte Carlo procedure used here, necessarily 
include approximations to make the models sufficiently tractable to a study of 
the type that is the subject of this report. Those approximations, although 
yielding appropriate macroscopic behavior for a nonequilibrium gas, do not 
adequately describe its microscopic behavior. For example, molecular 
encounters can occur that violate energy and momentum conservation. One, 
therefore, is at a loss as to how to treat the negative energies and momentum 
that arise during a simulation. 

Semiclassical models (refs. 10, 15, and 16) appear to be based on physi- 
cally realistic criteria; however, the Pearson and Hansen model was too sim- 
plified. Although this model would not seriously violate equilibrium concepts, 
it was subject to slow drifts from equipartition (ref. 10) (i.e., the rota- 
tional temperature would drift from equality with the translational tempera- 
ture) . Although, by heuristic arguments this model could be altered to 
qualitatively satisfy proper interaction behavior, we have based our investi- 
gation described here on the semiclassical model of Itikawa (ref. 16) . This 
model is founded on more rigorous concepts and, in addition, allows for treat- 
ment of molecular collisions. The model also satisfies conservation of 
probability. 
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In this paper we treat translation-rotation interactions for a uniformly 
flowing gas far removed from solid boundaries. In fact, we assume that the 
rotational relaxation does not affect the flow; we are concerned only with 
understanding the rotational relaxation behavior. Results are given based on 
calculations starting from three different sets of initial conditions: 

(1) complete equilibrium, (2) nonequilibrium equipartition (distribution of 
energy is constant for all molecules - the energy value is based on the total 
energy within the rotational states at complete equilibrium that is to be uni- 
formly distributed to each molecule; that is, total energy in rotation is 
physically correct but distributed incorrectly), and (3) nonequipartition 
(same as (2) except partition of energy within rotational states is not based 
on complete equilibrium) . Throughout this paper the key mathematical rela- 
tions essential to the Monte Carlo simulation are defined. Readers desiring a 
more comprehensive treatment of the method are referred elsewhere (refs. 1 and 
6-9) . 


PROCEDURE 


The essential aspects of the procedure are described in the introduction. 
The analytical relations peculiar to this investigation are described in the 
text that follows. 


Selection Rule Defining the Occurrence of an Encounter 

The key to an accurate simulation is the procedure for selecting the 
molecular pair to reach in a collision, determining whether a reaction occurs, 
advancing the time parameter in a systematic manner until the next collision 
occurs, and so on. Both the probability and time intervals are strongly 
dependent on collision frequency which, in turn, is dependent on the uncertain 
relation for intermolecular potential. 


To aid this discussion, it is worthwhile to refer briefly to classical 
relations and how those relations depend on, for example, an intermolecular 
potential base and on a two-parameter repulsion between the point-center model. 


In this case, where the intermolecular potential V(R) is spherically 
symmetric, the deflection angle of an encounter, X(b,g), which depend on 
impact parameter b and relative velocity of approach, g, and the &th moment 
"transport" cross section, qW, are given by 


X(b,g) 


! X 




(bdR/R 2 )//l - (b/R) 2 


- V(R)/(l/2)yg 2 


( 1 ) 


(g) = 2 t r f (1 - cos Z x)bdb 


( 2 ) 
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where Rq and y are the distance of closest approach and reduced mass, 
respectively. (For example, see chap. 8, ref. 12.) The collision frequency 
v is then given by 


V = nQWg (3) 

From this relation, we can compute the collision time At of an encounter and 
the elapsing time t. There results 


At 


2 1 
N v 


t = ^2 At 


(4) 


where N is the number of particles in a simulated cell. The collision 
frequency is not, however, accurately known in general. 

In the case when the potential can be described by 

V (R) = ~ (5) 

R 6 


the frequency is given by 

v 


If v max and Smax are the 
sionless ratios in the equation 

= M- 

v max/ \Smax 

define a curve for specified values of 6. If one accepts the above relation 
as representative of molecular encounters, then any point in the region below 
the curve represents a valid encounter, and points above the curve are invalid. 
We can then use the relation as a "selection rule," defining whether an 
encounter occurs or not. 

Rather than accept the above relations as completely valid, we also 
investigated results using the linear relation 


k (6-4)/6 


(7a) 


. 2/6 

6d \ „(£).., (5-4) / 6 

to | | A'' (S)g'- 


( 6 ) 


maximum values possible in a cell then the dimen- 


v 

v max 


B + (1 - B) — 

§max 


(7b) 
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where B is an adjustable parameter that gives the best results in the case 
B - 0.3. The reasons underlying this choice are described later. 

Discussion of transition probabilities is given in the next section. 
Given the fact that a rotational transition has occurred, however, trajecto- 
ries are required by the method so that the particle coordinates in phase- 
space can be recomputed. 


Collision Dynamics 

The relative velocity after the collision is obtained by knowing the 
rotational energy and momentum before and after a collision. These relations 
are classical relations given by 

(g 1 ) 2 = g 2 “ “ C^ri ~ E rl + E r2 - E r2 ) (^ 

and 

[gb - (Mr! - M rl + m; 2 - M r2 )/y] 

b’ = (9) 

g 

where E r2 and M r2 are rotational energy and momentum before a collision 
and a prime denotes value after a collision, the transition probabilities as 
well as calculation of a trajectory are based on |g ! - g|/g <<c 1 and 
|b* - b|/b << 1; that is, the relative velocities and impact parameters are 
only slightly perturbed as a result of the rotational transitions. 

The deflection angle (eq. (1)), although also dependent on intermolecular 
potential, can be adequately approximated to obtain trajectory results by the 
infinite-rise, rigid-spherical molecule of diameter a. Such a relation 
yields 

X(b) = 2 cos -1 (10) 

To properly account for the effect of inelastic collisions, we use the 
following mean value: 


X(b) - * 2 *< b ')J 


( 11 ) 


The derivation of the relation for rotational transition probabilities is 
briefly reviewed in the next section. 


ROTATIONAL TRANSITION PROBABILITY 


The details of the method and its applicability are discussed in 
reference 16. Only major derivations and results will be presented here for 
the purpose of direct application to the Monte Carlo simulation. 
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The interaction considered in the calculation of rotational transition 
probability is described by the reaction: 

N 2 C j i D + N 2 (j 2 ) ^N 2 (j|) + N 2 (j 2 ) (12) 

The calculation is based on a semi-classical theory. To make numerical calcu- 
lations tractable, several approximations are made. Whether these approxima- 
tions are valid is difficult to assess, except that they lead to the correct 
qualitative behavior in the several applications considered. The total wave 
function of the system is expanded in terms of a set of wave functions based 
on a "rotational" Hamiltonian and given by: 

^ =2 c a( t )l a> ex p(- rv) (13) 

a \ n / 

Here a specifies the rotational state of the molecules and W a and |a> are, 
respectively, the eigenenergy and eigenfunction of that state. The time- 
dependent coefficient C a (t) is then determined by: 

dC 

ift -gjir- = £<<*' l V l a > exp(iw a , a t)C a (14) 

where V is the interaction potential and u a i a = (W a , “ W a )/h. We separate 
the Hamiltonian into isotropic (spherically symmetric) and nonisotropic parts. 
Since our interest is in inelastic collisions, the elastic process due to the 
isotropic part of the potential, vt°), is eliminated by introducing the 
distorted-wave type coefficient defined by: 

D a (t) = C a (t) exp (| f l V (0) [R(t')]dt’} (15) 

where R(t) is the distance separating gravity centers of the molecules at 
time t. The coefficient D a satisfies the equation 

iK^=E[<a'|v|a> - v(°)6 a , a ]exp(im a , a t)D a 

= £<“' I V | cx) exp(ico a , a t)D a (16) 

a 

and the initial condition 

r l for a = a Q 

D a (t = -oo) = - (17) 

0 for a ^ a Q 

The second line in equation (16) defines the reduced matrix element (ot'jvla). 
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The transition probability for the process a Q -> a is given by 

P(a 0 •> a) = |D a (t = co)| 2 (18) 

In a rigorous treatment, the state a depends on the rotational 
angular momenta, j, and the projections, m, of both molecules. For the 
present problem, we are only interested in the probability averaged over the 
m states. We use Rabitz method (ref. 17) to eliminate the in-dependence of 
the interaction matrix (effective potential method) . We solve equation (16) 
with: 


a = D a 


D, i 
J 1^2 


and 


<06 ' | V | 0t> = ^ 



This treatment is discussed in more detail in appendix A. 

We can further approximate the solution of equation (16) by introducing 
the exponential approximation (refs. 18-20). 


D a' (“) = (<*' I exp K|a Q ) 

OO 

2 <a ' I K n I «o> 


n=o 


(19) 


where K is a matrix whose element is defined by 

* OO _ 

(ai|K|otj) = - i J dt^ajv^la..) exp (iaiQ^ t) (20) 

The element of K n is evaluated by 

<a'|K n |a 0 )= <a' |K|a n _ 1 )<a n _ 1 |K|a n _ 2 ). . .(axlKlao) (21) 

“1«2> • • • > a n -i 


For N 2 + we select the following interaction potential as the 
relavant interaction (refs. 21-22), 


V = v(°) (R) + (R) [P 2 (cos Xi) + P 2 (cos X 2 )] 

+V*- 2 -' (R)P 2 (cos Xi)P 2 (eos X 2 ) 


( 22 ) 


Here P 2 is the Legendre function of order two, R is the separation distance 
between the two molecules as before, and Xi is the angle between the 
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directions of the inter-molecular vector and the axis of ith molecule; the 
two molecules are not necessarily in a single plane, v(^) induces a simultane- 
ous rotational transition in both molecules in first order. Each term of the 
potential is assumed to have the form 


V (0) (R) = C exp(-ctR) - C 6 /R 6 ' 
V ^ (R) = AC exp(-aR) 

V^fR) = BC exp(-aR) 

i 

(See appendix A for further discussion.) 


(23) 


In the present calculation the trajectory R(t) is determined by solving 
the classical equation of motion approximated by 



-aR ^6 Eb 2 
C e + — r- - — o 

R C 6 Rc 2 


C24) 


where y is the reduced mass of the system, b is the impact parameter, Rq 
is the distance of closest approach, and E is the kinetic energy of the 
relative motion. This provides a semiclassi cal version of the modified wave 
number approximation (ref. 23) with the assumption that the gradients of the 
attractive part of the potential and the centrifugal force are much smaller 
than those of the repulsive potential, in the region where most of the rota- 
tional transition takes place. The solution of equation (24) can be obtained 
analytically and the time integration in equation (20) can be performed 
readily. We have, then 


( a i|R|“j) = -i ( a i I a j ) I 

with 


(25) 


and 


I 


TT 

olA 


cosech 



A = K(2mE)" 1/2 
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Two important properties of the transition probability should be 
mentioned: the conservation of probabilities and detailed balancing. The 

exponential approximation self-ensures the conservation of probabilities: 

X) p <Jij 2 + w) = 1 w 

w 

To satisfy the detailed balancing relation, we choose, as the kinetic energy 
E in the calculation of <a^ | Veff | aj>, the mean value of the initial and 
final channel energies (i.e., E = (1/2) (E^ + Ej)). The channel energy E^ for 
the ith channel is defined by 

Ei = E t - W a . = (E 0 + W aQ ) - W a . (27) 

where Eq is the initial relative kinetic energy and E t is the total energy 
(which is conserved during the collision). This procedure results in a 
symmetry relation 


p Uij 2 - j 1 , j 2 l ;E 0 ) = P(j 1 «j 2 t j 1 j 2 ;E 0 ? ) (28) 

where 

E 0 + Wj lj2 = E 0' + Wjj'jj,' 

In order to have a properly detailed balance, we modify our result by 

(2j , < + 1) (2j 9 < + 1) 

p (jli 2 + jl’j 2 ') = (2j 1 + 1) (2j 2 + If- PCjlj 2 " j i' j 2 ,) 

for (j^) ^ Cj 1 'j 2 '), and J (29) 

PCiii 2 ^ iij 2 ) = 1 - ^2 P(iii 2 ji , j 2 ') 

ji'j 2 ’ i j 2 ) 


Here j^< designates the smaller of and These modified transition 

probabilities, P, satisfy both the detailed balancing relation and the 
conservation of probabilities. 

The program listings for calculating rotational transition probability 
are presented in appendix B. 
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RESULTS AND DISCUSSION 


Calculations were performed using three different sets of initial condi- 
tions. These initial conditions, described in the introduction, are: (1) com- 

plete equilibrium, (2) nonequilibrium equipartition, and (3) non-equipartition. 
The first case, equilibrium, was run to test whether the method remains in 
equilibrium for long computational times; that is, to verify that the model 
would not drift to improper internal distribution (ref. 10) . The second case 
tests whether, indeed, the model has an internal mechanism to drive itself to 
equilibrium within a reasonable physical time scale. The third case provides 
insight into relative internal time scales to reach (1) a quasi -Boltzmann 
distribution characterized by a rotation temperature T rot T and then 
(2) the time scale for this quasi-distribution to decay to equilibrium 
Trot = T. The simulations, therefore, permit us to observe energy partition- 
ing and the relaxation mechanisms, as well as relaxation rates. Because we 
are here interested only in rotational transitions that lead to a final 
equilibrium state, all elastic collisions have been ignored in order to 
expedite the calculations. 


Equilibrium Case and Collision Frequency 

Many of the physical properties for one-dimensional calculations have 
been based on the hyperbolic function trajectory sech(at) (ref. 15); some 
others are based on the effective potential, and the classical equations of 
motion are solved for this interaction potential including the step by step 
energy conservation. 

Collision frequency, however, cannot be evaluated analytically for these 
potentials (see eqs. 1-3). It is also not feasible to use the Monte Carlo 
method to compute and thereby describe intermolecular potentials numerically. 
Nevertheless, we can semi- empirically determine a macroscopically (statisti- 
cally) correct collision frequency by using the results of equation (7) for 
the equilibrium case; that is, we determine the most probable index of power 
6 in equation (7) by varying 6=4 (Maxwellian molecules) to 6 = » (hard- 
sphere molecules) . The simulation should maintain equilibrium when the proper 
value of 6 is chosen. The effect of the parameter 6 on rotational energy 
is sensitive to the higher velocity collisions. We find that rotational tem- 
perature increases as the parameter 6 increases. Figure 1 shows the com- 
puted rotational temperature history from an initial equilibrium state as the 
value of 6 is varied. The proper value for which the model (point centers 
of repulsion) fits close to the present model, is found to be approximately 
C = 0.45, or 6 =7.3, where C = (6 - 4)/6. Also presented in the same fig- 
ure are the results of the temperature variation obtained by using the linear 
selection formula of equation (7b) where the parameter B was varied from 0 
to 0.5. The best value of B that satisfies this selection rule seems to be 
B = 0.3. Temperature histories computed using this value (C = 0.45) are shown 
next . 


In figures 2 through 7, translational and rotational temperatures for 
different initial conditions are plotted as functions of nondimensional 
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collision time. These temperatures should asymptotically approach the 
equilibrium values shown by the dotted line at the end of the time scale. 

Also shown in figure 2 for several time intervals, t/T = 0, 5, 10, and 20, are 
the translational and rotational distribution functions. The translational 
distribution function ftr is plotted as a func tion of molecular speed ratio 
c/c Q , where c Q is the rms speed defined by /3kT/M. The rotational distribu- 
tion function f r is presented in terms of the rotational energy level j 
for the same collision times corresponding to the translational distribution 
functions. Mean collision time t is calculated based on the initial trans- 
lational temperature, kept always at room temperature Ttr - 320° K. Note 
that initial equilibrium distribution functions for rotational and transla- 
tional energies are selected at random from the Rayleigh and Maxwellian dis- 
tributions, respectively, (their distributions, therefore, do not represent 
analytical functions), and all subsequent distributions evolve from the 
present simulation using the Itikawa model. 

Throughout the testing, all distributions and temperatures represented 
equally valid equilibrium states. Figures 3a and 3b show the comparison of 
the Monte Carlo solutions with theoretical (equilibrium) functions for the non- 
dimensional time at t/x = 15 and 20. Since the Monte Carlo solution - the 
so called Klimontovich function (ref. 24) - does not, in general, present a 
smooth function, a mean distribution function (Boltzmann solution), time 
averaged over the last five collision times (from t/x = 15 to 20), is shown 
in figure 3c. This mean value can be compared with the theoretical value when 
the gas approaches equilibrium. The agreement with theory seems to be very 
good. 


Nonequilibrium Case 

The result of the equipartition , nonequilibrium test is shown in 
figure 4. The gas initially starts in the equipartition state with the rota- 
tional energy of all molecules assigned the mean equilibrium energy correspond- 
ing to the level of j = 10, and where the translational energy is specified 
as in the previous case. Several energy distribution functions are shown at 
nondimensional times of 0, 1, 2, 3, 5, 20, and 35. The simulation again seems 
very good. The distribution functions for t/x = 30 and 35, and the mean dis- 
tribution functions averaged over these time intervals, are compared with 
theoretical calculations in figure 5; once again the agreement is good. 

Notice that only even-numbered rotational energy levels are occupied. This 
follows since we have considered a homonuclear model to be initially in an 
even-numbered level, and we have not allowed changes in nuclear spin. 

The last test investigated is for the nonequipartition and nonequilibrium 
case. All gases initially start with the rotational energy level of j = 12 
(corresponding to a rotational temperature of T r = 455° K) and where the 
translational energy is selected at random from the Maxwellian distribution 
corresponding to a translational temperature of 320° K. In figure 6, both the 
rotational and translational temperatures approach an equilibrium state, 
corresponding to a temperature of 374° K. This occurs in an exponential 
manner. Transition distribution functions are shown at t/x = 0, 1, 2, 3, 5, 
and 20 in the same figure. The present model also appears to perform very 
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well in this case. This is the case (initial rotational temperature higher 
than translational temperature) for which the modified Pearson-Hansen model 
failed to show satisfactory performance when extended to higher temperatures 
(ref. 10) . 

The distribution functions for t/T = 15 and 20, and the mean distribu- 
tion functions of rotational and translational energy, averaged over t/x = 15 
to 20, are shown in figure 7. Comparison with theoretical distributions is 
very satisfactory. 

Notice that in both figures 4 and 6 the gas relaxes asymptotically to 
equilibrium. The approach is rapid initially and quite slow finally. 

Of all the models tested to date, the semiclassical model of Itikawa 
appears most satisfactory. However, this model with 1000 molecules consumes 
about 100 sec of the CDC 7600 machine computing time to travel one character- 
istic collision time. The model, therefore, may still need further simplifi- 
cation to permit its practical use in the more complex molecular simulations, 
such as, for example, shock wave structure and gas-surface interactions. 


CONCLUSION 


The Itikawa model when used with the appropriate representations for 
collision frequency provides an adequate physical description of a homonuclear 
diatomic gas in rotational relaxation. This model appears to hold the most 
immediate promise for further application to more complex problems. However, 
because of the computational time required to do rigorous calculations at each 
Monte. Carlo collision event, subsequent investigations will require simplifi- 
cation of the algorithm. Nevertheless, the present model will permit one to 
examine the principal features of rotational effects in nonequilibrium flow of 
diatomic molecules, such as shock wave structure. 


Ames Research Center 

National Aeronautics and Space Administration 
Moffett Field, Calif., 94035, July 1975 
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APPENDIX A 


INTERMOLECULAR POTENTIAL AND INTERACTION MATRIX ELEMENT 
Interaction Potential 

Each term of the interaction potential (eq. (23)) is assumed to be 
defined by the representations given in the text. The first term, V™), can 
be determined fairly reliably from either the result of molecular beam experi- 
ments (ref. 21) or the analysis of transport coefficients (ref. o 22). These 
parameters are C = 3.44xl0 3 eV, Cg = 73.4 eV A -6 and a = 3.16 A -1 . Unfortu- 
nately, we have little information about the anisotropy of the interaction 
potential for N 2 + N 2 . Therefore, we adopt the form given by vC 1 ) and v( 2 ) 
in equation (23), similarly as in reference 15, and regard A and B as adjust- 
able parameters. The values employed in this report are A = B = 0.2. 


Interaction Matrix 

Applying the effective potential method (ref. 17) to the interaction 
potential given by equations (22) and (23), we can calculate the matrix 
element as follows: 

<j 1 'j 2 '|V eff |j 1 j 2 > = <j 1 'j 2 '|V eff |j 1 j 2 > - V (0) (R)6 jl , jl 6 ; j 2 , j2 

= CJ 1 * |ji3 2 ) c exp [-aR(t) ] 


(Al) 


where 


(V 

V o 


H 

o 


+ 1) C 2 j 2 + l)(2j 1 ' + 

1) (2 j 2 ' + 1)] 1/4 C-1) L 


[(-l) ;i2 (2j 2 + 1)' 1/2 

( h o SK'J, 


jl (25 1 + ir 1/2 ( j2 ’ 

0 o) S jl'n] 


j 1 ' jl 2\/3 2' 

0 0 oA o 

o J) E C^l) 1,2 (o 

m=o ,2,4 

2 ^ 
0 0/J 



(A2) 


^ is the 3 - j symbol and L = max (j 1 +j 2 j j j ' + j 2 ') 
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APPENDIX B 


PROGRAM LISTING FOR ROTATIONAL TRANSITION PROBABILITY CALCULATION 


*DFCK MAIN 

PROGRAM MAIN! I NPUT, OUT PUT , TAPF5 =T NPUT ,T APE6=0UT PUT ) 

CROSS SECTION CALCULATION 

COMMON /LLMIT/ LL 1M I N , L L1M AX , LL 2M I N ,Ll 2MAX 
DIMENSION PWAVE<20,20) , SI GMA { 20 , 20 ) ,LABCSA( 16) 

DIMENSION SIG(20) 

COMMON /CM1/ PWAVE ,E KI N , LI 0 , L20 ,B I MP, NMAX , NPR I NT, LiPAR, L2PAR 

COMMON /CM2/ J EL 1 , J EL 2 , L ABCSA , L LMAX 

COMMON /CM3/ VA , VB , I PR T1 , I PRT2 ., I PRT 3 , 1 PRT 4 

COMMON / C MVR1 / VC , V ALPHA , VC 6, BB IMP, EEEE 

COMMON /CV1/ ABCC (AO ,40,9 ) 

JMAX= 20 
P AI2= 6.283185 

C** INITIAL CLEAR OF ABCC, SHOULD BE MADE FOR EACH VA,VB 
DO 5009 1 =1,40 
DO 5009 J= 1 » 40 
DO 5009 K=1 , 9 
ABCCU,J»K)= 0.0 

5009 CONTINUE 

5010 CONTINUE 

DO 5016 I = 1 , J MA X 
DO 5016 J=l, JMAX 
S I GMA ( I , J ) =0. 0 
PWAVE( I » J ) = 0.0 
5016 CONTINUE 

C** TNPUT ** IMPACT PARAMETER (IN ANGSTROM) 

C** CALCULATION FROM 'BIMPI' TO ' B I MP F ' WITH STEP * D B I M P * 

READ( 5 » 4 ) BIMPI ,OBIMP, BIMPF 
4 F ORMAT ( 3F 10.0) 

IF(BIMPI.LT.O.O) GO TO 5999 

C** TNPUT ** CROSS SECTION CALCULATION (ICROS=i) OR NOT (0) 

RE AD( 5, 3 ) ICROS 

C** INPUT ** RELATIVE K IN FT IC ENERGY (IN EV) 

READ (5,1) EKIN 

1 FORMATt F10.0) 

C** TNPUT ** INITIAL ROTATIONAL STATES 
READ( 5,2) LiO , L20 

2 FORMAT (215) 
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C** INPUT ** MAX NO. OF TERMS IN EXP 
R E AO (5 , 3 ) NMAX 
3 FORMAT U 5} 

C** INPUT ** INDEX FOP PRINT OUT 

RE AD (5, 6 ) IPRTI, IPRT2, IPRT3, IPRT4 

6 FORMAT (415) 

C** INPUT ** POTENTIAL PARAMETERS FUR SPHERICAL PART 
C** V ( R J = VC*EXP (-VALPHA*R)-VC6/R**6 

C ** V IN EV, R IN ANGSTROM 

R E AD (3,7) VC, VALPHA, VC6 

7 FORM AT (3F10.0) 

C~ A INPUT ** POTENTIAL PARAMETERS FOR NGN-SPHERICAL PART 
READ (3,5) VA , VB 
5 FORMAT ( 2 F5.0 ) 

C** INPUT ** WHEN IPKJSG = 0, PRINT PARTIAL SUM OVER B I MP 
R EAD( 3,3) IPRTSG 

C ** INPUT ** LIMITATION OF RANGE OF L i , L 2 
C*-* L I = L 1M I N-L L MAX , L2=L2MIN~L2MAX 

READ (5,8) LLMIN, LIMA X,l 2M I N , L2MAX 
3 FORMAT ( A I 5 ) 

L L 1MI N= L1MIN+I 
L L 1M AX= L1MAX + I 
L L 2 M I N = L2MIN+1 
L L2MAX = L2MAX+I 
J l M I N = (LLIMIN + D/2 
J 1 MAX= (LLIMAX+D/2 
J2MIN= (LL2MIN+U/2 
J2MAX= ( LL2MAX+1 ) /2 
WRITE (6, 5 900) 

1 RIMPI tORIMP,BIMPF, 

2 I CRUS , 

3 FKIN, 

L 1 0 , L 2 0 , 

5 NMAX, 

ft IPRTi, IPRTi» t IPRT3, IPR T4, 

7 VC , VALPHA ,VC6 , 

3 VA , VB , 

9 IPRTSG 

9900 FORMAT! 1H1/5X » 10HINPUT DATA/ 

L 3X, L8H8I MPI , DR IMP , fci I MP F / , 3 F 10 .2/ 

2 5X» 6HICR0S/ , 15/ 

3 5 X , 5HEK I N/ , F 1 0 .6 / 

4 5X , 8HL L0,L20/,2I 5/ 

5 3 X , 5 HNMAX / ,15/ 

6 *3 X , 24H IPRTI , I P RT2 » I PRT3 , 1 PRT4/ , 4 f 5/ 
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7 5X,i4HVC,VALPHA,VC6/,3F10.4/ 

8 5X » 6HVA ,VB/,2F5. 2/ 

9 5X, 7HIPRTSG/, 15// ) 

WRITEI6 ,5902) L1MIN, UMAX, L2MIN, L2MAX 
5902 FORM AT (1H0/5X,4H L1=I4,3H - ,I4,3X,4H L2=I4,3H - ,14/) 
£*•*«=<<**** NPRINT 

NPR I NT = NMAX-3 
B IMP* BIMPI 

IFUCROS.EQ.O ) IPRT 4 = 0 
IFUCROS.EQ.O) GO TO 5500 
H3= DB IMP/3. 0*PAI2 
NSUM* 0 
MPRINT=0 
H3B* H3* B I MP 
CALL PROS 
5 LOO CONTINUE 

00 5111 I * JIM IN , J IMA X 
00 5111 J=J2MIN, J2MAX 

S I GMA( It J> = S I GMA ( 1,0) +H3B *P WA V£ ( I , J) 

5111 CONTINUE 

B IMP=BIMP+DBI MP 
CALL PROB 
H4B*4.0*H3*BIMP 
00 5112 I = J1MI N» J1MAX 
DO 5112 J=J2MIN, J2MAX 

S IGMA( I »J)=SIGMA( I , J ) +H4B*PW AV F ( I,J ) 

5112 CONTINUE 

BIMP* BIMP+DBIMP 

CALL PROB 

H3B*H3* B I MP 

DO 5113 I = JIM I N , J IMA X 

00 5113 J* J2MIN, J2MAX 

S IGMAC I , J) = SIGMA! I , J ) +H3 B* PWAV E ( I , J) 

5113 CONTINUE 

IFUPRTSG.5Q.1) GO TO 5120 
GO TO 5300 
5120 CONTINUE 

I F ( BIMP. GT.BIMPF-0.5*0BIMP ) GO TO 5200 
N SUM* NSUM+2 
GO TO 5100 
5200 CONTINUE 
MPRINT *1 

IFIIPRTSG .EQ.O) GO TO 5400 
5300 CCNTINUE 

C***************SfcARCH FOR MAX OF SIGMA 
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SMAX=0.0 

DO 5305 I=J1MIN, JlMAX 
DO 5305 J= J2MIN, J2MAX 

IFtSIGMAI I,J) .GT.SNAX) SMAX = SIGMAII,J) 

5305 CONTINUE 
r************ PRINT 

L 1 PARO= MODI L10 ,2 ) +1 
L2PAR0= MODI L 20, 2 ) + 1 

WRITE 1 6 ,5901 ) BIMPI , BI MP, NSUM , L 10, L 20 , EKIN 
5901 FORMAT! 1H1/5X ,23HCR0SS SFCTION FOR BIMP=F6.2,3H - ,F6.2,7H (NSUM= 
1I2,1H),3X, 5HI L 10= I 2, 1X,4HL 20=1 2, 1H) ,3X,5HFKIN=F 13 . 5/ /3X ,2HL2 ) 

00 5311 JINV=J2MIN, J2MAX 
J= J2MIN + J2MAX-JINV 

L 2= 2*J-3+L2PAR0 

C ********** NORMALIZATION TO SMAX 
DO 5310 I = JIM I N , J IMA X 
SIGtI)= S IGMAI I, JI/SMAX 

5310 CONTINUE 

JF= J1MIN + 15 

IFIJF.GT .JIM AX) J F= JlMAX 

WRITE! o, 5911) L2, I SIG! I) , I=J1MIN, JF ) 

591 L FORMAT! 1H0,3X, 12, 16! 1PF6.3, IX) ) 

5311 CONTINUE 

WRITE! 6, 5912) I LABCSA! I ) ,1=1 ,16) 

5912 FORMAT! 1H0,4X, 16! 3X, 12, 2X) , 3H LI) 

S MAX= 0.1*SMAX 
WRITE! 6, 5919) SMAX 

5910 FORMAT ! 1 H0/3X, 15HN0RMALIZ FACTOR, 2X ,~ 13 . 5 , 2 X, 1 1HA NGS TROM* * 2) 

SEL= -SIGMA! JEL1,JF:L2) 

C********* ELASTIC 

WRITE(6,5918) SEL 

5918 FORMAT ! 1 HO »4X ,2 9HT OT AL INELASTIC CROSS S c CT ION , 2X , E 1 3 . 5 ) 

1 F (MPR I NT . EQ. 1) GO TO 5400 
GO TO 5120 

5400 CONTINUE 

C********* ROTATIONAL CONSTANT IN EV 
BRCT= 0.2512E-3 
S = 0.0 
SWW= 0.0 

W0= BROT * FL0AT(L10*!L10+1)+L20*!L2 0+1) ) 

DO 5415 I=J1MIN, JlMAX 
L 1= 2*1— 3+L1P AR 0 
DO 5414 J= J2M IN , J2MAX 

IFtl.FQ. JEL1.AND. J.EQ. JFL2) GO TO 5414 
L 2= 2*J-3+L2PAR0 
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W 1= BROT * FLOAT (L1 J ML1+1) + L2*(L2+11 ) 

W= <W0-W1)**2 
SHW= Srt W + S IGMA ( I , J 1 *WW 
S=S+SIGMA ( I ,J) 

541 A CONTINUE 
5415 CONTINUE 

WRITEI6, 59161 S 

5916 FORMAT! 1HO , 4X , 2 9HTOTAL INELASTIC CROSS SECT ION, 2X , E13.5 1 
WRITE(6, 5917 1 SWW 

5917 FORMAT ( 1H0/4X»27H** 3 $ £ **' SUM OF PCORCT ** ***, 3x ,4HSW W = E 13 . 5/ 1 
GO TO 5010 

5500 CONTINUE 

IF(BIMP.LT.O.O) GO TO 5010 
CALL PROB 

IF(BIMP.GE.8IMPF-0.5*DBIMP)G0 TO 5010 

5501 CONTINUE 
BIMP=BIMP+QBIMP 
GC TO 5500 

5999 CONTINUE 
STOP 
END 


*DECK PROB 

SUBROUTINE PROB 

C********* main PROGRAM FOR THE CALC. OF TRANS. PROB. 

COMMON / L L M I T / l L 1M I N , L L 1M AX , LL 2MI N , LL 2MA X 
01 HENS I ON PODD ! 20,20 1 , PEVN (20,20 ) ,P WAVE <2 0, 20 ) 

0 I MENS I ON AKSUMO! 20 , 20) , A K SUM1 ( 20 ,2 0) 

0 I HENS ION LABCS A ( 16 ) 

COMMON / M VI/ AMATRX(20,20,9) 

COMMON /MV2/ VBB ,VAA ,RROT, ETOT ,BRC , VVALP 
COMMON / CMVR1/ VC , V ALP HA , V C6 , B IMP , E EEE 
COMMON /MV3/ NCOUNT 

COMMON / C Ml/ PWAVE ,E KIN , H0,L20 ,BBB B, NMAX ,NPRINT, L1PAR,L2 PAR 
COMMON / CM2/ J5L 1 , JE L2 , LABCS A, L LMAX 
COMMON /CM3/ VA , VB , I PR T 1 , I PRT2 , I PRT 3, I PRT4 
I C LOCK= 0 
B I MP=BBBB 
J MAX=20 

LMAX= 2* JMAX-2 

r********* REDUCED MASS IN AMU 
RMASS= 14.02 

C ********* ROTATIONAL CONSTANT IN FV 
BROT = 0. 2512E-3 


18 



VVALP= 0 . 04-5723# V AL P HA/ SURT(RMASS) 

VAA = 0.4472136*VA 
V BB= 0.2*VB*0. 6298283 

ET0T= EKIN + FLOAT (L10#(L10+1)+L20*(L20+1) )#RROT 
EB = ETOT/BROT 
E Bl= SQRT(EB) 

LLMAX= I NT {E B 1 ) +1 

IF { IPRT1.EQ. I) GO TO 11 

W RITE! 6 , 901) B I MP , L 1 0, L20 , EKIN , ETOT ,LLMAX,EB, VC ,VALPH A , VC 6 , VA , VB 
901 FORMAT (1H1 »5X»5HBIMP=F6.2» 5 H ANGST » 5 X »6 H ( L 10= 1 2 , 5H L 20= 12 » 2H )» 

1 5X,5HEKIN=F 13. 5,3HEV ,3X,5HETOT=E 13.5 »3X,6HLLMAX= 12, 

2 3X,3HEB=F8.2/6X» 3HVC= F 10 • 2 » 2X» 7HVALPHA=F 8. 3, 2X ,4HVC6=F 7. 2 , 2X , 

3 3HVA=F7.2»2Xf3HVB=F7.2/> 

11 CONTINUE 

I F (LLMAX. GT.LMAX ) L L M AX = L MAX 
IF(LLIMAX.GT.LLMAX) LL1 MAX = LLMAX 
IFILL2MAX.GT. LLMAX) LI. 2MAX= LLMAX 
DC 18 I = 1 1 JMA X 
DO 18 J = 1 » JM A X 
P EVN ( I » J ) = 0.0 
PODDII t J) = 0.0 
PWAVEII , J) = 0.0 
AKSUMO ( I , J)=0.0 
AKSUM1 ( I ♦ J ) = 0.0 
DO 16 M= 1*9 
A MATRX ( I , J,M)= 0.0 
16 CONTI NUF 
18 CONTINUF 

L L 10= L 10 + 1 
L L20= L20 + 1 
J F L 1= (LL10+D/2 
JEL2=(LL20+l)/2 
L 1PAR= MOD ( LI 0 » 2 ) + 1 
L 2PAR= MOD ( L20» 2 ) + 1 
I FIL1PAR.LT. LUMIN) L1PAR= LL1MIN 
IF(L2PAR.LT.LL2MIN) L2PAR= LL2MIN 
C L0= FLOAT ( ( 2*L 10+1 )*( 2*L20+ l) ) 

PE0= 0.0 

C ********* n- i 

N = 1 

L LI I = LL10 
L L2I = LL20 
NCCUNT = 0 
DO 59 K= 1 » 9 

CALL VMATRXILL LI »LL2I »K) 
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ii in ■ 


59 CONTINUE 

LL1FMN= MAXOILL 10-2.L1PAR) 

LL1FMX= M I NO ( LL 10+2, I.L 1MAX ) 

LL2FMN= M AXO < LL20-2 » L2 PAR ) 

L L2FMX= MIN0ILL20+2, LL 2MAX ) 

K = 1 

JJU= (LL1I + D/2 
J J I 2= (LL2I + l)/2 

00 64 1 = 1,3 
LL2F= LL2 0-4+2* I 
DO 63 J= 1 , 3 
LL1F= LL 10—4+2* J 

IF(LL2F.LT«1. 0R.LL1F.LT .1 ) GO TO 62 
J JF1= (LL1F+1J/2 
J JF2= ( LL2F+1 )/2 

AKSUMOl JJF1, JJF2) = A MATRX ( J J 1 1 , JJ 1 2 «K ) 

62 K =K + 1 

63 CONTINUE 

64 CONTINUE 

T CT PW= 0.0 
T CTPC= 0.0 

DO 74 LL1F=LL1FMN,LL1FMX,2 
DO 73 LL2 F=LL2FMN» LL2FMX,2 
J JF 1= (LL1F+D/2 
JJF2= ( L L 2F + 1 ) / 2 

PODDIJ JF 1 , JJF2)= AKSUMOIJ JF1,JJF2) 

PEVNIJJFl , JJF2)= 0.0 

1 FUL1F.EQ.LL II .AND . LL 2F .E Q .LL 2 1 ) P EVN ( JJ F 1 , J JF 2) = 1 . 0 
PWAVE( J JF1,JJF2)= PODD ( JJ FI , J J F2 ) ** 2*4 .O+PE VN ( J JF 1 , J JF 2 i **2 
P=PWAVF( JJF1,JJF2) 

LO 1=LL 1 F- 1 
L02=LL2F - 1 

I F ( IPRT 2.F.Q. 1 ) GO TO 72 
PCORCT = P 

WRITE (6, 902) N, L01 , L02 , P, P CORCT , PODDi J JF1 , J JF 2 ) , AKSUMOl J J F 1 , J JF 2) 
902 FORMAT! 1H0,2X,2HN=I 2, 3X,5H L1F = I2,2X,5H L2F = I 2 , 5X ,6HPW AVF = E13 .5 , 

L 3X,7HPC0RCT=E13.5,3X,5HPnDD=Ei3.5, 5X,6HAKSUM=E 13.5) 

72 CONTINUE 
T0TPW=T0TPW+P 

73 CONTINUE 

74 CONTINUE 

IF(IPRTl.EQ.l) GO TO 75 

WRITE(6» 904) N, TOTPW, ICLOCK.NCOUNT ,TOTPC 

75 CONTINUE 

I F(NMAX.EQ.l) GO TO 300 
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N = 2 

C******$**ioo 

100 CONTINUE 
NC CUNT * 0 
N2= 2*N- 2 

LL2 1 MN=MAXO ( LL20-N2 » L2 P AR } 

LL2IMX= LL20+N2 

IF(LL2IMX.GE.LLMAX) LL2IMX= LLMAX 
DO 149 LL2I=LL2IMN»LL2 !MX»2 
J J 1 2= (LL2I+1I/2 

E B2= SQRTIEB- FLOAT ( LL 2 I- l ) **2 ) 

LLMAX1= 1 NT ( E B2 ) +1 
L L II =LL 1 0+N2 
JJli=(LLlI+i)/2 

IF(LLMAXI.GE.LLIMAX) L LMAX1= LL LM AX 
I F(LL1I .GT .LLMAXU GO TO 130 
DC 114 K=5 t tit 3 
CALL VMATRXILL1I ,LL2I ,K) 

114 CONTINUE 

DC 119 K = 3 »9 t 3 

CALL VMATRXILLlt, LL2I ( K) 

119 CONTINUE 

IFILL1I-2.LT.1) GO TO 125 

A MA TR X( JJI ltJJI 2t4) = AMA TR X < J J I 1 - 1 , J J 1 2 t 6 ) 

K =4 

I F(LL2I .F0.LL20-N2.0R.LL2 I . EQ .LL20 +N2 ) CALL VM ATRXI L L 1 1 , LL 2 I t K ) 
AMATRX( J JI 1, JJI 2t 7 ) = AMATR X < J J 1 1- 1 , J J I 2 + 1 , 3 ) 

K = 7 

I F (LL2I . EU.LL20+N2-2 .0R.LL2 I ,EQ.LL20*N2 ) CALL VMATRXt LL II , LL 2 I ,K) 
IPILL2I-2.LT. 1) GO TO 130 

A MATRX (JJI It JJ 12, I )= AMATR X< JJI 1- i , J J I 2-1 , 9 ) 

K = 1 

IFILL2I . FQ.LL20-N2+2.0R.LL2I .EQ.LL20-N2) CALL V MAT RX ( L L 1 1 , LL2 I » K J 
125 CONTINUE 

IF(LL2I-2.LT.l) GO TO 130 

AMATRXl JJI l,JJl2t2) = AMA TR X( J J 1 1 , J JI 2- 1 1 8 ) 

K= 2 

I F (LL2I .EQ.LL20-N2 ) CALL VMATRX ( L L 1 I t LL2 I t K) 

130 CONTINUE 

L L 1 ! -LL 1 0-N2 
JJU* (LL1I+D/2 
I FILL1I.LT. L1PAR) GO TO 149 
DO 134 K = 5 ,8 t 3 
CALL VMATRXILL1I ,LL2I ,K) 

134 CONTINUE 
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DO 139 K= 1*7,3 

CALL VMATRX(LL1I,LL2I ,K) 

139 CONTINUE 

AMATRXt JJI 1,JJI2,6J= AMATRX ( JJ 1 1 +1 , JJ 12 ,4) 

K = 6 

IF(LL2I.EQ.LL20-N2.0R.LL2I.EQ.LL20+N2) CALL VMATRXt LL1 I , LL 2 I ,K) 
AMATRX (JJIl»JJI2»9)= A MATRX ( J J 1 1+ 1, JJ 12+1,1 ) 

K = 9 

I F(LL2I .EQ.LL20+N2-2.0R.LL2I .PQ.LL20+N2) CALL V MATRX ( LL 1 1 , LL2 I , K) 
IF(LL2I-2.LT.1> GO TO 149 

AMATRX( JJI l ,JJI2 ,2) = AMATRX (JJI 1 , JJ 12-1 ,8 ) 

K= 2 

IFILL2I.EQ.LL20-N2) CALL VMATRX ( LL 1 1, LL 2 1 ,K ) 

AMATRX! JJI 1 , J JI 2 , 3) = AMATRX (JJI 1+1 , JJI2-1 ,7) 

K = 3 

I F(LL2I. EQ.LL20-N2+2.0R.LL2I . C Q .LL 20-N2 ) CALL V MATRX( LL II , LL 21 , K) 
149 CONTINUE 
C********* 150 

LL1 1 MN= MAX0(LL10-N2+2,L1PAR) 

L L 1 1 MX= LL10+N2-2 

I F { LL1 IMX.GE.LL1MAX) LL1IMX= LL IMA X 

DC 199 LL 1 I = LL1JMN»LL1 I MX, 2 

JJI1 = (LL1I + U/2 

EB2 = SORT (EB- FLO AT ( LL 1 1- 1 )**2 ) 

LLMAX1=INT (EB2 ) +1 
L L 21 = LL20+N2 
J J 12= (LL 2 I + 1 )/2 

I F (LLMAX1 .GE.LL2MAX) L LMAX 1= LL2MAX 
IFJLL 2 I.GT .LLMAX1) GO TO 180 
DO 164 K =5 , 9 

IFtLLll .EU.LL10+N2-2.AND.K.FQ.6) GO TO 163 
CALL VMATPX(LL1 I,LL2I,K) 

GO TO 164 

163 CONTINUE 

AMATRX ( JJI 1 , JJ 12, 6 )= AMATR X ( JJ 1 1+ 1 , J J 1 2 ,4 ) 

164 CONTINUE 

IF (LL2I-2.LT . 1) GO TO 175 

AMATRX(JJI1»JJI2»2)= AMATR X ( J J I 1 , J J 1 2- 1 , 8 ) 

AMATRXt JJI 1,JJI2,3)= A MATRX ( J J I 1+ 1 , J J 1 2"1 , 7 ) 

I F (LLl 1-2 »LT . 1 ) GO TO 180 

AMATRXt JJ II, JJ 12,1 )= AMATRX (JJI 1-1, JJ 12-1-9) 

175 CONTINUE 

IFtLHI — 2.LT.1) GO TO 180 

AMATRXt JJI 1,JJ 12,41= AMATRX ( J J 1 1- 1 , JJ 12 , 6 ) 

180 CONTINUF 
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LL2I=LL20-N2 
J J 1 2= (LL2I + 1J/2 
IFILL2I .LT.L2PAR) GO TO 199 
DO 184 K = 5,6 

IFILLiI.EC).LLiO+N2-2.AND.K.EQ.6) GO TO 183 
CALL VMAT RX (LL1I»LL2I»K) 

GO TO 184 

183 CONTINUE 

AMATRXl JJI1, JJI2,6)= AMATRX l JJ I 1+ 1, JJI2,4) 

184 CONTINUE 

DO 189 K= 1 ,3 

CALL VMATRX I L LI I *LL2I»K) 

189 CONTINUE 

AMATRXl JJU,JJ 12,8 ) = AMATRXIJJI l, JJ I 2+ 1,2) 
AMATRXl JJIL,JJI2,9)= AMATRX ( JJ I 1+1, JJ 12 +1 , 1 ) 
IF(LLlI-2.Lr. 1) GO TO 199 

AMATRXl JJI i,JJ 12, 4) = AMATR X ( J J 1 1- 1 , JJ 1 2 , 6 ) 
AMATRXI JJI1,JJI2,7)~ AMATRX I J J 1 1- 1 , J J 1 2 + 1 , 3 ) 

£********* igg 

199 CONTINUE 
LL1FMX=LL10+2*N 
LL2FMX= LL20+2*N 
L L1FMN= LL10— 2*N 
IFILLiFMN.LT. L1PAR) GO TO 201 
LL 11=LL 1FMN+2 

200 CONTINUE 
LL2FMN = LL20-2*N 
IFILL2FMN.LT. L2PAR) GO TO 202 
L L22= LL2FMN+2 

GO TO 205 

201 CONTINUE 
LL1FMN= L 1 PAR 
L L 1 1= L 1 PAR 
GO TO 200 

202 CONTINUE 

L L 2FMN=L2 PAP. 

L L22= L2PAR 

205 CONTINUE 
LL22F= LL2FMX-2 

IF! LL2FMX • LE . LL2MAX ) GO TO 206 
LL2FMX= LL2MAX 
L L 2 2F = LL2MAX 

206 CONTINUE 

DO 249 L2=LL22,LL22F ,2 
L L 1 1F= LL1FMX-2 



EB2® SQRTIEB- F LOAT < L2-1 ) **2 ) 

L LMAX1= I NT ( EB2 ) +1 

IF(LLMAXl.GE.LLiMAX) LLMAXl= LLiMAX 

IF(LLlFMX.GE.LLMAXi) LL11F= LLMAX1 

00 248 LI= LL11,LL11F,2 

Ki* 1 

K 2= 0 

L22=L2-2 

Li 1= Li-2 

IFIL2.LF.L2PAR) L22® L2PAR 
I F(Ll.LE.LiPAR) Lll® L iPAR 
IFIL2.LE.L2PAR) Ki=4 
IF(Li.lE.LlPAR) K2= 1 
K= Ki 

Jl® <Ll+l)/2 
J 2® (L2+D/2 
AK® AKSUMO ( Jl • J2 ) 

L2P2*L2+2 

DO 219 LL2® L22.L 2P2,2 
K= K+K2 
L lP2=Li+2 

00 218 LL1=L11.L1P2.2 
JJi® (LL1 + D/2 
JJ2® (LL2+1J/2 

AK$UMi{ J Jl, JJ2) = AK*AMATRX( J1,J2,K) +AKSUMK JJI , JJ2) 

K= K+l 

218 CONTINUE 
2 19 CONTINUE 

c********* PRINT of AM AT R X 

IF (N.LT. I PRT3) GU TO 248 

L01=L1-1 

L02® L2- 1 

IF (L2.GE.LL22F-1.AND. I L20 + 2*N.LE. LLMAX+2) GO To 220 
I F ( Li »GE .LL liF-1 • AND • L L10+2*N. LE. LL MAX + 2) GO TO 220 
IF(L2.FQ.LL22.AND.LL20-2*N.GE.L2PAR-2) GO TO 2 20 
IF{Ll.EQ.LLll.AND.LL10-2*N.G£.LiPAR-2) GO TO 220 
GC TO 248 
220 CONTINUE 

WR ITE I 6 » 906 ) L01,L02 

906 F0RMAT(1H0,3X ,15HCHECK OF AMATRX » 5X »3HL1=I2,2X» 3ML2® I 2 1 
WRITE (6, 907) (AMATRXI J1»J2, J) , J=1 ,9) 

907 FORMAT ( 1H ,2X,9E13.5) 

248 CONTINUE 

2 49 CONTINUE 

0********* 250 
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LL1FMX= LL10+2+N 
T OTPW= 0.0 

00 299 LL2F=LL2FMN,LL2FMX,2 
L L 1 1 F= LL1FMX 

E 82= SQRT(EB- FLOAT (LL2F-1)**2) 

LLMAXi= I NT( E82 J +1 

I F ( LLMAX 1 . GE.LL1MAX) LLMAX1= LL IMA X 
IF(LHFMX.GE.LLMAX1 ) LL1 1 F= LLMAX1 
DO 298 LL1F=LL IFMN,LL11F,2 
J JF 1= (LL1F + D/2 
JJF2 = (LL 2F+l)/2 

AKSUMM J JF1,JJF2) = AKSUM1 ( JJF1 , JJF2 )/ FLOAT (N) 

I F ( MOD (N»2).NE.O) GO TO 294 

AKSUMK J JF1,JJF2) = -4.0*AKSUM1 ( JJFl ,JJF2) 

PEVNUJFl, JJF2> = PEVNl JJF l, JJF2 )•«- A KSUM1 ( J JF 1 , J JF2 ) 

PWAVE( J JF1 , JJF2 ) = PODD ( JJF It J J F2 ) **2*4. O+PE VN (JJF1,JJF2)**2 
GO TO 295 

294 CONTINUE 

PODDUJFl , JJF2>= PODD( JJFl, JJF2)+ AKSUM 1 ( J J F 1 , J JF2 ) 

PWAVEl J JF 1 .JJF 2) = POOD (JJFl , JJF 2 )** 2*4.0+ PEVNt JJFl , JJF2)**2 

295 CONTINUE 

P = PWAVEt JJF1,JJF2) 

L 0 1= LL IF - 1 
L02= LL2F-1 

C********* PRINT OF PWAVE 

I F (N.LE . NPRI NT ) GO TO 297 
I F (P.GE.0.1E-4) GO TO 296 
IF(LLIF.EQ.LLIFMN) GO TO 296 
I F(LLLF.GT .LL11F-2) GO TO 296 
GC TO 297 

296 CONTINUE 
PCORCT = P 

W R 1 TE ( 6 » 902) N , LO 1 , L02 , P, PCORCT , POD D ( J JF1 ,JJF2) ,A KSUM1 ( J JF 1, J JF2) 

297 CONTINUE 
TQTPW= TOTPW+P 

AKSUMOt J JF l, J JF2) = AK.SUM1 ( JJFl , JJF2 ) 

AKSUMK J JFL, JJF2 )= 0.0 

298 CONTINUE 

299 CONTINUE 

PE i=PW AV E ( JEL I, JEL2 > 

PE 10= ABS( (PE1-PE0) /PEI ) 

1 F ( I PRT 1 . EQ. 1 J GO TO 3904 

WRITE (6, 904) NtTOTPW , I CLOC K , NC OUNT , TOTPC 
904 F0RMAT(1H0//2X,8H***** N=I 2 » 5X ♦ 13HSUM OF PW AV E= F.l 3 . 5 , 
15X,7HICLQCK=I 10»3HSEC,5X, 7HNC0UNT = I 10, 5X, 6HT0TPC=E13. 5 ) 
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3904 CONTINUE 

I F1ABS1T0TPW-1.0) .LT.0.1E-3.AND.PE10.LT.0.1E-3) GO TO 300 
I FtN.EQ. NMAX) GO TO 300 
N=N+ 1 
P E0= PEI 
GC TO 100 

300 CONTINUE 

C ;********* ELASTIC 

DELST1= PEVN! JEL1, JFL2) 

DELST2* -2.0*P0D0( JFL1, JEL2) 
c ********* FINAL PRINT 
MPR I NT = 0 

L 1PAR0= M0D(L10,2H-1 

301 CONTINUE 

WRITE! 6, 901) BJ MP, L1C, L20 ,EKIN, ETOT ,LLMAX ,EB, VC ,V ALPHA , VC 6 , V A , V B 
WRITE! 6, 910) N,TOTPW,I CLOCK, TOTPC 

910 FORMAT! 1H ,5X ,5HNMAX= 12 ,2X ,6HT0TPW= Eli .5, 2X .6HCL0CK* 15, 3HSEC, 

1 2X,6HT0TPC=E13.5//3X,2HL2) 

I MN= LL2FMN 
I MX=LL2F MX 

I F !MOD! LL 2FMX-L2PAR , 2) .NF . 0) I MX* LL2FMX-1 
DO 319 1= IMN, IMX, 2 
L L2F= IMX+IMN-I 
J 2= (LL2F + 1J/2 
LL11F* LL1FMX 

EB2= SQRT(EB- FLOAT! LL2F-1 )**2) 

LLMAX1* I NT ( EB2 ) +1 

I F ! LLMAX 1 . GE . LL 1MAX ) LLMAX1* LL1MAX 
I FtLLlFMX.GE. LLMAX1) L L11F =LLMAX1 
J 1 MX* (LL11F + U/2 
J IPX 1= J 1 MX 
J1MN= (L1PAR+D/2 

I F! J1MX-J 1MN.GT . 15) JiMX* J IMN ♦15 
L02=LL2F-1 

WRITE(6» 911) L02,(PWAVE!J1,J2) , J 1 = J IMN , J1 MX ) 

911 FORMAT! 1 HO/ 3X , 12, 16 ( 1PF6.3 , IX) ) 

I F(MPP INT.EQ. 1) GO TO 319 
LMIN2* M I NO! L 20 , L 02 ) 

CLIO* FL0AT(2*LMI N2+1) 

DO 318 Jl= J 1 MN , J IMX 1 
L01= 2*Jl-3+L 1PAR0 
L MINI* MI NU( L10 , L01 ) 

C L 1= CLIO* FL OA T ( 2* LM I N 1 + 1 ) 

P= PWAV E 1 J 1, J2 ) 

PCCRCT* CL1/C LO*P 
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TOTPC= T QT PC+PCORCT 
PWAVE ( J 1 » J2) = PCORCT 

318 CONTINUE 

319 CONTINUE 

L ABCSA ( 1 ) = LI PAR — 1 
DO 320 1=2,16 

L ABCSA ( I)=LABCSA(1)+(I-1)* * * * 2 

320 CONTINUE 

W R I TE( 6, 912J (LABCSA! J) ,J=1,16) 

912 FORM AT (1H0»4X»16(3X» 12, 2X) ,3H LI) 

IFIMPRINT.EQ. 1) GO TO 9906 
c********* PRINT OF PCORCT WITH PELASTIC MODIFIED 
PWAVE! JF.Li ,JEL2) = 1 . O-TOTPC + PW AVE < J EL 1 , JE L 2 ) 

IF! IPRT4.EQ. 1) GO TO 1000 
M PR INT= 1 
GO TO 301 
9906 CONTINUE 

WRITE16, 9909 ) TOTPC 

9909 FORMAT! 1H0/2X , 35H*******PC CRCT WITH ELASTIC MODIFIED, 5X , 

1 7H( TOTPC =E13.5, 2H )/) 

1000 CONTINUE 

PWAVE! JE LI ,JEL2) = PW AVE ( JE L 1 , JEL2 ) - 1 .0 
WRITE!6, 9911 ) DEL ST 1 ,DE LS T 2 

9911 FORMAT! 1H0, IX, 28H***** FOR ELASTIC SCATTER ING , 3X, 3HD 1=E 13 . 5, 
l 2X, 3HD2=E 13. 5) 

RETURN 

END 


*DECK V MAT 

SUBROUTINE VMATRX! LL1 , LL2 , K ) 

COMMON / M VI/ AMATRX!20,20,9) 

COMMON /MV2/ VBB , VAA ,BROT , ETOT , BRC , VVA LP 
COMMON / C MVR1 / VC, V ALPHA, VC6 , B I MP , E IJ 
COMMON / M V3/ NCOUNT 
COMMON /CVl/ ABCC (40, AO, 9) 

JJ1= (LL1 + D/ 2 
JJ2 = (LL2+1I/2 
L 1 J= LLi-1 
L 2 J= LL2-1 

I F ( K.GT . 3 ) GO TO 1011 
1.21= L2J-2 

IFIL2I.LT. 0) GO TO 1510 
GO TO 1100 

1011 IF (K.GT. 6) GO TO 1012 
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L 21 = L2J 
GC TO 1100 
1012 L2 1 = L2J+2 
1100 CONTINUE 

IF(M0D(K,3).NE.i) GO TO 1111 
L 1 1 = L1J-2 

IF(LII.LT.O) GO TO 1510 
GO TO 1200 

1111 I F ( MOD (K»3).NE.2) GO TO 1112 
L 11= L 1 J 
GO TO 1200 
11 12 L 11= L 1J + 2 
1200 CONTINUE 

ABC= ABC C. ( LL1»LL2»K) 

IF (ABC .NE. 0.0) GO TO 1299 
CALCULATION OF V EF F ( L 1 I » L 2 I / L1J,L2J ) 

CC1= CG20 ( Lil ,L1 J) 

C C 2= C G 2 0 ( L 2 1 ,L2J) 

1211 CONTINUE 

I F ( K.Lt • 5 ) CS IGN= ( -1 . 0 ) ** ( L 1 J + L2 J ) 

IF(K.GE.6) C S IGN= ( - 1. 0 ) ** ( L 1 1 + L2 1 ) 

C = FLOAT ( (2*L1I + 1)*(2*L2I + 1)*( 2*L 1J +1 ) * ( 2*L 2J + 1 )) **0. 25 *CSIGN 

1212 CONTINUE 

B= VBB*CC1*CC2 

1213 CONTINUE 
A = 0.0 

I F(L2I.NE.L2J ) GO TO 1215 
A= VAA*CC1/ SQRT ( FLOAT (2*L2I+1)) 

I F ( MOD (L2I*2).NE.O) A= -A 
1215 IF(LII.NE.LIJ) GO TO 1219 

A A= VAA*CC 2/ SORT (FLOAT (2*L1 1 + 1 ) ) 

I F ( MOD (LlI f 2).NE.O) AA=-AA 
A = A+AA 
1219 CONTINUE 

ABC= C*(B+A) 

ABCC(LL1 ,LL2»K)= ABC 
1299 CONTINUE 

W I = BROT * FLOAT ( L 1 1 *( L 1 1 +1 ) +1.2 1 * ( L 2 1 +1 ) ) 

WJ = BROT* FLOAT (LI J*( HJ + 1 ) + L2J*(L2J + i> ) 

W I J= ABS( WI-WJ) 

E 1= ETOT-WI 
EJ= ETOT-WJ 

IF(EI.LE.O.O.OR.FJ.LE.O.O) GO TO 1510 
EIJ= 0.5*(EI+EJ ) 

CALL ROOT ( RC ) 
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BRO 1 .0-(BIMP/RC )**2+VC6/EIJ/RC**6 
E I J1 = E I J* BRC 
I F ( K .EQ • 5 ) GO TO 1500 

IF(L1I.EQ.L2J • ANO.K . EQ .3 ) GO TO 1500 
I F(L1I.EQ.L2J. AND.K. EQ.7) GO TO 1500 
A I J= VVA L P* SQRT (EIJD/WIJ 
6W = EIJ1/WIJ 
D B ALPH= EW/AIJ 
A P A I = 1. 57C796327/A I J 
F = EXP(-APAI) 

F A I J= 2.0*APAI*F/(1.0-F*F) 

A AA= ABC#DBALPH*FAIJ 
C*** 999 CHECK PRINT 
959 CONTINUE 

1498 CONTINUE 

NCOU NT = NCOUNT + l 

1499 CONTINUE 

AMATRXt J J 1 , JJ2,K )= AAA 
RETURN 

1500 CONTINUE 

A AA= A8C/VV41P* SQRT(EIJl) 

GO TO 1498 
1510 CONTINUE- 
A A A= 0.0 
GO TO 1499 
END 


SUBROUTINE ROOT(RC) 

C***** REVISED 8/26/74 

COMMON / C MVR1 / VC , V ALP HA , V C 6 , B IMP , E Kl N 
R M I N= 1.12 

RRO= ALOG( VC/EKIN) /VALPHA 
IF (RR0.GE.4. 1) R R 0= 4.1 
N = 1 

3099 CONTINUE 
RR=RRO 

3100 CONTINUE 

RL= EKIN*BIMP**2/PR**2 

V = VC* EXP(-VALPHA*RP) 

V 1 = — V AL P HA*V 
IFIVC6.EQ.0.0) GO TO 3101 
VR= VC 6/ RR**6 

V = V-VR 

V 1= V1+6.0*VR/RR 
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3101 CONTINUE 

F = (V+RL-EKlN)/ (2.0*RL/RR-V1) 

IFIABS(F/RR).LT.0.1F-5) GO TO 3199 
I FIN. GE. 100) GO TO 3299 
R R = RR+F 

IFIRR.LT. RMIN) GO TO 3900 
N=N+1 

GO TO 3100 
3199 CONTINUE 
RC= RR 
RETURN 

3299 CONTINUE 

WRITE (6» 998) RR,F,RRO 

998 FORMAT (1H0//5X, 14HERR0R N GT 100 , 3X ,3HRR= F 1 3. 5 » 3X f 2E 13. 5/ / ) 
RR=RRO 
GO TO 3199 
3900 CONTINUE 

RR0= 0. 5*< RRO+RMI N) 

GO TO 3099 
ENO 


*DECK CG20 

FUNCTION CG20 ( J 1 » J2 J 

I FIJ2.EQ.JI+2 .QR.J2.FQ .Jl-2) GO TO 8001 
IFIJ2.EQ.Jl) GO TO 8002 
C = 0.0 
GO TO 8100 

8001 CONTINUE 

IF ( J2.EQ. Jl+2) J = J1 
IFIJ2.EQ.Ji-2) J= J2 
X 1= FLOAT (J«-2)/FL0AT (2*J+5 ) 

X 2- FLOAT ( J+ 1) /FLOAT { 2*J+ 3 ) 

X3= l.O/FLOAT (2*J+1) 

C = SQRT I 1.5*X1*X2*X3) 

GO TO 8099 

8002 CONTINUE 
J= J1 

X 1= FLOAT! J+i) /FL0AT(2*J+3) 

X2= FLOAT ( J ) /FLOAT ( 2 *J +1 ) 

X3= 1.0/FL0AT(2*J-l ) 

C= - SQRT (X1*X2*X3) 

8099 I F ( MOD (J»2).NE.G) C= -C 

8100 CG20= C 
RETURN 
END 
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Figure 1.- Variation of rotational temperature depending on specified values for the potential 

parameters. 





Relaxation time, t/r 
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Figure 2.- Temperature and distribution functions: Complete equilibrium initial conditions. 





(a) t/x = 15 

Figure 3.- Translational and rotational distribution functions: Complete 

equilibrium initial conditions. 
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